function [Optima] = Niching(RegionSet,SampleSet,RegionSampleId,Radial)
%Nitching: obtain the final optima set
% RegionSet: ~ -by- (2+~) matrix, the whole subregion set
%         [partition depth, partitionable, features ...]
% SampleSet: ~ -by- (1+d) matrix, the original sample set
%         [objective value, x1, x2, ...]
% RegionSampleId: ~ -by- ~ cell, the id of the samples in each region
% Radial: the radial distance that be regarded as neighbor

%% the possible optimal solutions
N = size(SampleSet,1); % number of samples
OptCandidate=zeros(N,1);
for i=1:length(RegionSampleId)
    if RegionSet(i,2)==0
        OptCandidate(RegionSampleId{i})=1;
    end
end
%% ranking the solutions
[~, Ranking] = sort(SampleSet(:,1));
SampleSet = SampleSet(Ranking,:); % sort the samples from bad to good
OptCandidate = OptCandidate(Ranking,:);
%% clearing
OptimaSet = OptCandidate; % belongs to the optima set or not
sampleid = 1;
while 1
    Distance = (sqrt(sum((repmat(SampleSet(sampleid,2:end),[N,1])-SampleSet(:,2:end)).^2,2)))';
    neighbor = find(Distance<=Radial);
    [~,best_neighbor] = min(SampleSet(neighbor,1));
    OptimaSet(neighbor) = 0; % delete all neighbor
    OptimaSet(neighbor(best_neighbor)) = 1; % save the best neighbor to the optima set
    if sampleid == (neighbor(best_neighbor))
        next = find(OptimaSet((sampleid+1):end),1);
        if isempty(next)
            break
        else
            sampleid = sampleid + next;
        end
    else
        sampleid = (neighbor(best_neighbor));
    end
end
Optima = SampleSet(OptimaSet==1,:);
end
